arXiv: 1507.04375vl [astro-ph.SR] 15 Jul 2015 


Preprint typeset using style emulateapj v. 05/12/14 


A UNIFIED COMPUTATIONAL MODEL FOR SOLAR AND STELLAR FLARES 

Joel C. Allred 

NASA/Goddard Space Flight Center, Code 671, Greenbelt, MD 20771 


Adam F. Kowalski 

Department of Astronomy, University of Maryland College Park 
AND 

Mats Carlsson 

Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029, Blindern N-0315, Oslo, Norway 


ABSTRACT 

We present a unified computational framework which can be used to describe impulsive flares on 
the Sun and on dMe stars. The models assume that the flare impulsive phase is caused by a beam 
of charged particles that is accelerated in the corona and propagates downward depositing energy 
and momentum along the way. This rapidly heats the lower stellar atmosphere causing it to ex¬ 
plosively expand and dramatically brighten. Our models consist of flux tubes that extend from the 
sub-photosphere into the corona. We simulate how flare-accelerated charged particles propagate down 
one-dimensional flux tubes and heat the stellar atmosphere using the Fokker-Planck kinetic theory. 
Detailed radiative transfer is included so that model predictions can be directly compared with ob¬ 
servations. The flux of flare-accelerated particles drives return currents which additionally heat the 
stellar atmosphere. These effects are also included in our models. We examine the impact of the 
flare-accelerated particle beams on model solar and dMe stellar atmospheres and perform parameter 
studies varying the injected particle energy spectra. We find the atmospheric response is strongly 
dependent on the accelerated particle cutoff energy and spectral index. 
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1. INTRODUCTION 

Stellar flares are the result of large explosions in the 
atmospheres of stars. They are produced when mag¬ 
netic helds, which have been stressed by convective mo¬ 
tion in the stellar photosphere, reconnect rapidly releas¬ 
ing their stored energy. In addition to directly heating 
the reconnection site in the corona, much of this en¬ 
ergy goes into accelerating charged particles to high ve¬ 
locity. These travel down magnetic field lines colliding 
with the increasingly dense plasma, depositing their en¬ 
ergy and momentum, and quickly heating the plasma in 
the reconnecting flux tubes to temperatures > 20 MK. 
The high temperature causes emission to dramatically 
brighten in virtually all regions of the electromagnetic 
spectrum. Flares are ubiquitous and have been observed 
on stars of nearly all spectral types. 

Because of the Sun’s close proximity, understanding 
solar flares is especially important. Together with coro¬ 
nal mass ejections, these directly affect the Earth envi¬ 
ronment. They have significant impact on spaced-based 
communications, the power grid, and the manned space 
program. Also because of the Sun’s proximity, it is possi¬ 
ble to spatially resolve many of the features which drive 
solar flares. Such observations are not possible on other 
stars. However, despite the lack of spatial resolution, 
there is still much to be learned by studying stellar flares. 
(In this paper, we will refer to solar flares as those ex¬ 
clusively on the Sun. Stellar flares will refer to flares on 
all stars except the Sun.) Active dMe stars are known 
to have flares much larger than those on the Sun. These 
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stars Spin faster and generate much larger magnetic fields 
(jHawlev et al.l[2?)l^ and references therein) resulting in 
flares some 10 ^ times more energetic than typical so¬ 
lar flyes le.g.. iH awley fc PettersenllI99H iHawlev et al.l 
1200, H iKowalski et al.ll2010( l. In addition, the background 
intensity of dMe stars is much smaller than that of the 
Sun. Thus, when the se stars flare, the signa l can be much 
higher. For example. IKowalski et al.l (1201011 observed the 
optical luminosity of the dMe star, YZ CMi, to brighten 
by ~ 200 times. In comparison, in the largest solar flare 
an increase in irrad i ance of just ^100 ppm was observed 
(j Woods et al.ll2005 iMoore et al.ll20l4l . 

Accelerated electrons are known to play an important 
role in transporting energy during flares. Their presence 
can be detected from the bremsstrahlung radiation they 
produce as they collide with the ambient plasma. Since 
these electrons are a major source of flare heating, it 
is crucial that we accurately model them. Fortunately, 
the Rama ty High Energ y Solar Spectroscopic Imager 
(RHESSI: iLin etffl 1200211 observes this bremsstrahlung 
radiation from which it i s possible to deduce the injected 
electron spectrum (e.g.. iHolman et all 1200311 . Thus, to 
simulate how the lower atmosphere responds to this heat¬ 
ing, we model the precipitation of these electrons from 
the acceleration site in the corona to the footpoints in 
the chromosphere and below. 

In addition to electrons, ions a re also likely accelera ted 
by reconnecting magnetic fields. lEmslie et al.l (|20I2ll es¬ 
timated the energy in accelerated ions to be comparable 
to that of accelerated electrons in many flares. How¬ 
ever, since the bremsstrahlung cross-section is inversely 
proportional to the square of the mass of the collid- 
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ing particle (|Haug|[r997f) . their presence is much harder 
to directly detect. Even with this limitation, RHESSI 
has detected th e presence of ycelerated ions in severa l 
large flares le.g.. lHurford et aD2nnfiHEmslie et al.ll2ni2[l . 
Even though direct evidence of their presence is scarce, 
models of particle acceleration predict that lower energy 
ions (< 1 MeV) will b e accelerated -albeit on longer time 
scales than electrons (|Petrosian fc Liul[2004ll . Therefore, 
to accurately model the flaring atmosphere the effects of 
flare-accelerated ions must also be included. 

A beam of charged particles propagating down a mag¬ 
netic tub e induces an electric fi eld which drives a return 
curre nt (Ivan den OordI 119901 : IZharkova &: Gordovskvvl 
I2006L and references therein). The magnitude of the cur¬ 
rent depends upon the particle beam spectrum. In fact, 
the return current alters the particle spect rum causing 
a flattening at low energies (|Holmanl l2012fl . This cur¬ 
rent additionally heats the ambient plasma through Joule 
heating. In the flaring corona, this can be a major source 
of energy and could explain superhot coronal tempera- 
tures (> 30 MK) observed in sever al large flares (e.g. 
iCasrii fc Linll20iri iCasni 6111111201411 . 

It is the purpose of this paper to present a unified com¬ 
putational model which can describe the atmosphere of 
both solar and stellar flares including the processes most 
important to flare dynamics. To do this we simulate the 
transport of a beam of non-thermal particles injected at 
the top of a magnetic flux tube and follow the subse¬ 
quent heating of the stellar atmosphere. Our model flux 
tube extends from the sub-photosphere into the corona. 
Pressure, temperature, and density vary by many orders 
of magnitude across these regions, and our models must 
be able to accurately represent these very different con¬ 
ditions. For example, in the corona radiative transfer is 
dominated by numerous high temperature, non-LTE, op¬ 
tically thin atomic transitions. In the photosphere, how¬ 
ever, radiative transfer is optically-thick and in LTE. In 
the chromosphere, neither of those approximations hold 
and the full radiative transfer equation must be solved 
in detail for several important atomic transitions. Flares 
produce high-speed shocks (i.e., > 600 km s“^). We 
have found that to resolve these requires a computational 
grid with spacing < 100 m. Currently, models which in¬ 
clude detailed radiative transfer at such high-resolution 
are only tractable in one-dimension. 

Decades of observations have revealed that flares cer¬ 
tainly have a three-dimensional geometry. However, the 
strong magnetic forces present in flaring active regions 
confine charged particles to flow along field lines. Thus, 
to a good approximation flare dynamics can be mod¬ 
eled using a one-dimensional geometry with that dimen¬ 
sion being the axis of a model flux tube. To partially 
account for the three-dimensional nature of particular 
flares, the emission from individual flux tube models 
can be combined using timing and spatial information 
provided by observations (e.g., images from the Atmo¬ 
spheric Imaging Assembly on the Solar Dynamics Ob¬ 
servatory (SDO/AIA) and RHESSI). In this way, impor¬ 
tant processes are resolved in ways that are currently not 
tractable in fully 3D simulations. However, this method 
does not account for the interaction between plasma on 
differing flux tubes. This could affect heating rates since 
those depend on the plasma density which can be sen¬ 
sitive to mixing of plasma on differing tubes. We spec¬ 


ulate that the effect of mixing on plasma density will 
be small compared to that produced by chromospheric 
evaporation, which is captured by our ID simulations. 
This is because neighboring flux tubes are likely to have 
been heated by similar fluxes of non-thermal particles, so 
will have similarly elevated densities. In this paper, we 
present results from parameter studies varying ID model 
flux tubes and injected non-thermal spectra. In subse¬ 
quent papers, we will combine these ID models to form 
a more complete representation of particular flares. 

In Section [21 we describe our method to solve the 
equations of radiation hydrodynamics and the model¬ 
ing framework that we use to simulate solar and stellar 
flares. We discuss how thermal conduction and radiative 
transfer are included in the models. The chromospheric 
radiative transfer is of particular importance, and we 
describe our method to model emission from numerous 
optically-thick, non-LTE atomic transitions which dom¬ 
inate that region. We also describe how radiative back- 
warming from coronally produced X-rays and extreme 
ultraviolet (XEUV) radiation is included. In Section |3l 
we present how our models simulate the precipitation of 
flare-accelerated particle beams. We present our method 
for modeling the direct collisional excitation and ioniza¬ 
tion of the ambient plasma by these particles in Sectionj?) 
In the region of beam impact, these dominate over the 
thermal rates, significantly altering the radiative trans¬ 
fer. Section [5] describes our method for simulating how 
return currents additionally heat flaring flux tubes. In 
Section |6l we present results of a parameter study which 
we have conducted to understand the range in dynam¬ 
ics predicted from various injected particle beams. In 
Section |7l we summarize our results and present conclu¬ 
sions. 


2. RADIATION HYDRODYNAMICS 

Flares produce high-speed shocks and increase density 
and radiation throughout the stellar atmosphere. To 
model them, the radiative transport equation must be 
coupled with the standard equations of hydrodynamics. 
During flares, chromospheric plasma evaporates into the 
corona and waves can propagate into the photosphere. 
Modeling the atmospheric response to flare heating re¬ 
quires a model that can extend from the sub-photosphere 
through the chromosphere, transition region and into 
the corona. The transition region is extremely narrow, 
and we have found that accurately resolving it requires 
spatial scales as small as ^ 100 m. A model that in¬ 
cludes all of these elements represents a major compu¬ 
tational undertaking. We use th e RADYN co de devel¬ 
oped by iCarlsson fc SteinI (I1992L I1995L Il997f l to solve 
the equations of charge and population conservation cou¬ 
pled to the equations of radiation hydrodynamics in ID. 
We briefly summarize the method of solution here. The 
equations of radiation hydrodynamics are. 


dp dpv 
dt dz 


( 1 ) 
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dz 
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and p are 

the height, density, 

inter- 


nal energy density, velocity and pressure, respectively. 
g is the gravitational acceleration, and is a viscous 
stress term added to aid in achieving numerical stability. 
Fc and are the conductive and radiative fluxes, re¬ 
spectively. The conductive flux has the classical Spitzer 
form but is li mited so that it does not exceed the satu¬ 
ration limit of iSmith fc Aueil (|1980ll . In the atomic level 
population equation, rii is the number density in a given 
atomic state, and N' is the total number of atomic states 
considered in these simulations. We have derived a pa¬ 
rameter, mg, which represents the average mass of the 
plasma per hydrogen atom. This parameter assumes con- 
stant abundan c e rati os using the abundances derived by 
lAsplund et al.l (|2009ll . It has a value of 2.26 x 10“^^ g. 
Thus, p = monH, where uh is the hydrogen number den¬ 
sity. Pij is the transition rate from state i to state j and 
is given by Pij = Cij + Rij , where Cij and Rij are the 
collisional and radiative rates, re spectively. Th ese are 
discussed in much more detail in iMihalasI (|1978L Chap¬ 
ter 5). In the radiative transfer equation, is the 
frequency (v) and angle (/i) dependent specific inten¬ 
sity, and and Xvfj. are the emission and absorption 
coefficients, respectively. The radiative flux divergence, 
dFr/dz, is obtained by integrating Eq. [S]over frequency 
and angle. Qcor is a coronal heating term which is nec¬ 
essary to maintain a hot corona. Qbeam and Abeam are 
terms describing flare-accelerated particle beam heating 
and momentum deposition, respectively. They are de¬ 
scribed in more detail in Section |3l Qrc is a heating term 
due to return currents and is described in Section [SJ 
These coupled non-linear equations are solved implic¬ 
itly using a Newton-Raphson iteration scheme. Ad- 
vected quantities a re treated using the second-order up¬ 
wind technique of Ivan Fe ed (1197 711. RADYN employs 
an adaptive grid (|Dorfi fc Drurvlll98^ which is designed 
to resolve shocks and steep gradients which can form 
in flaring stellar atmospheres. The grid cell concentra¬ 
tion i s chosen to be propo rtional to the desired resolu¬ 
tion. iDorfi fc DrnrvI (jlQSTT l define a resolution operator 
which is strongly dependent on the absolute value of the 
gradients of the radiative hydrodynamic variables. Be¬ 
cause gradients of these variables can be large at different 
heights, weighting parameters have been chosen to give 
preference to those variables that require the most grid 
sensitivity. In these flare simulations, the largest weights 
were required on temperature, velocity, and atomic level 
populations, to properly resolve the transition region, 
strong shocks that form in the loops, and non-LTE pop¬ 
ulation densities which affect the convergence of the ra¬ 
diative transfer equation. Due to the complexity of the 
physical problem, we parameterize certain terms in our 
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Figure 1. The optically-thin radiative loss function used in these 
simulations. 

model using previous work. Thus, some uncertainties in 
our results may be only due to the uncertainties in these 
external models. 

2.1. Optically-Thick Radiative Transfer 

The RADYN code solves the radiative transfer equa¬ 
tion for the non-LTE conditions that dominate in the 
chromosphere. RADYN solves the atomic level popula¬ 
tion equations (Eq. 2]) for a six-level with continuum hy¬ 
drogen atom, a nine-level with continuum helium atom, 
a six-level with continuum Ca ii ion, and a four-level 
with continuum Mg ii ion. This allows the calculation of 
numerous transitions which are important to the chro¬ 
mospheric energy balance. Table [T] lists the line tran¬ 
sitions and Table [5] lists the continuum transitions that 
we model in detail. Eor these transitions, the radiative 
transfer equation is solved for up to 100 frequency points 
and five angular points providing us with detailed line 
profiles. 


2.2. Optically-Thin Radiative Transfer 

The densities in the transition-region and corona are 
typically low enough that the “coronal approximation” 
applies. In this case, the radiative transfer is dominated 
by numerous optically thin lines. These lines are formed 
when ions are collisionally excited and radiatively de- 
excited. Since the atmosphere is optically-thin in these 
regions, this radiation is assumed to escape, thus having 
a net cooling effect on the corona. We model the radia¬ 
tive transfer in these regions using a radiative loss func¬ 
tion which is obtaine d by summing all transitions in the 
CHIANTI database (iDere et al.lll997t iLandi et al.ll2013H 
except those already accounted for in Tables [U and [2j To 
generate this function the CHIANTI calculations were 
performed with the assumption of a constant electron 
density of 10^° cm“^. Eigure[T] shows this radiative loss 
function. 


2.3. XEUV Backwarming 

Half of the optically-thin radiative losses described 
above are directed outward and leave the stellar atmo¬ 
sphere. The other half are directed downward and will 
get absorbed in deeper and denser regions. This ad¬ 
ditional source of heating and ionization in the lower 
atmosphere becomes especially important during flares 
when the coronal XEUV emission c an become elevated 
by orders of magnitude. Previously, lAllred et al.l (120051 
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Table 1 

Bound-Bound Transitions 


Atom 

Xij (A)" 

Transition 

Atom 

Xij (A) 

Transition 

H I 

1215.67 

Lyo 

Ca II 

8662.14 

3d -H- 4p 2po^2 


1025.73 

Lyd 


8498.02 

3d ^D^/2 <-!■ 4p 


972.52 

Lyq 


8542.09 

3d ^ 1 ) 5/2 4p 


949.74 

Ly5 

He I 

625.56 

Is^ ^50 -H- Is 2s ^51 


6562.79 

Ho 


601.42 

Is^ ^ So 0 Is 2s ^ So 


4861.35 

H/3 


10830.29 

Is 2s ^Si e-l Is 2p 3pp 


4340.47 

Hy 


584.33 

ls2 1^0 0 Is 2p ipp 


18751.3 

Pao: 


20581.29 

Is 2s i5o Is 2p ipO 


12818.1 

Pa/3 

He II 

303.79 

Is ^5 'i/2 ^S-ij 2 


40522.8 

Bra 


303.78 

Is ^Si /2 0 2p 2pp 

Ca II 

3968.47 

H 

Mg II 

2802.70 

h 


3933.66 

K 


2795.53 

k 


^ These are vacuum (air) wavelengths for Xij below (above) 2000 A. 


Table 2 

Bound-Free Transitions 


Atom 

^ic 

Initial State 

Atom 

Me ('^) 

Initial State 

H I 

911 

n = 1 

He I 

504 

ls2 ISo 


3646 

n = 2 


2600 

Is 2s 


8204 

n = 3 


3121 

Is 2s ^50 


14584 

n = 4 


3421 

Is 2p ^P? 


22787 

n = 5 


3679 

Is 2p ipf 

Ca II 

1044 

4s 

He II 

228 

^‘S' 1/2 


1218 

3d ^ 733/2 


911 

2s ^‘S' 1/2 


1219 

3d ^P>5/2 


911 



1417 

4p 2pp^2 

Mg II 

824 

2p® 3s ^5 i/2 


1422 

"^3/2 


1168 

2p6 4p 2pp^2 




1169 

2p6 4p 2pp^2 


hereafter, AOS) developed a method to model how heat 
is deposited from XEUV backwarming, but that method 
did not account for the increased photoionizations from 
this flux. Here we present a new technique which self- 
consistently includes heating and photoionizations. We 
have used the CHIANTI database to tabulate emissiv- 
ities for numerous transitions as a function of temper¬ 
ature and wavelength. RADYN calculates the XEUV 
spectrum produced from a model loop by integrating 
the product of these emissivities with the emission mea¬ 
sures from the transition region and coronal portions 
of the loop. Finally, the radiative transfer equation is 
solved for several optically-thick chromospheric lines, as 
described above, assuming that this emission is incident 
on the chromosphere from above. In solving the radia¬ 
tive transfer equation, photoionization cross-sections for 
XEUV emiss ion are calculated and added t o the rates as 
described by iWahlstrom fc CarlssonI (1199411 . 

To understand how XEUV backwarming affects our 
loop models, we have generated loops which include 
and exclude this heating term. Figure [5] compares the 
structure of the QS.SL.HT loop model (see Section \m 
which has been generated with and without XEUV back¬ 
warming and using the technique of AOS. We find that 
the XEUV backwarming term results in a chromosphere 
which is 1000 - 2000 K hotter than would otherwise be 
expected. The technique presented here and that of AOS 
produce similar temperature structures. However, the 
radiative transfer predicted by these methods is signif¬ 


icantly different. This is illustrated in Figure |3] which 
compares the emission from the Ca ii H and He 1 10830 A 
lines for loops generated with and without XEUV back¬ 
warming and using the technique of AOS. Our technique 
produces a He i 10830 A line with a much deeper absorp¬ 
tion profile. This line is formed when continuum photons 
from the photosphere are absorbed by neutral helium 
atoms in the 2s ^Si excited state. The XEUV flux in¬ 
creases the photoionization rate of heliu m. These ions 
quick ly recombine into the 2s state (|Golding et al.l 
l2014f) . resulting in more absorption than would be ex¬ 
pected without the XEUV flux. The Ca ii can be un¬ 
derstood similarly. Without the inclusion of the XEUV 
flux, the region of the chromosphere where the Ca ii H 
line center forms is cooler resulting in an overabundance 
of ground state ions and an overall absorption profile. 
For both XEUV backwarming techniques, the Ca ii H 
line has a central reversal peaking in the near wings. 
The technique of AOS predicts more flux in the central 
reversal than our technique. Since their technique does 
not include direct photoionizations by the XEUV flux, 
it predicts an overabundance of Ca ii relative to Ca ill 
ions, and hence an increased flux in the Ca ii H central 
reversal. 

2.4. Opacity 

Of course, there are many continuum transitions not 
included in Tabled The contribution to the opacity due 
to these transitions is treated using the opacity package 
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Column Mass (g cm"’) 

Figure 2. The temperature as a function of column mass for loop 
models generated using the XEUV backwarming method described 
here (black), using the XEUV technique of A05 (green) and without 
XEUV backwarming (red). 


Ca IIH 



He I 10830 A 



Figure 3. Profiles for the Ca ii H and He i 10830 A lines from 
the loop models generated using XEUV backwarming (black), the 
technique of A05 (green), and without XEUV backwarming (red). 

of lGustafssonl This package constructs opacity as 

a function of temperature, density and frequency assum¬ 
ing that these transitions are in LTE. These are included 
as a background opacity source in the detailed calcula¬ 
tions of the transitions listed in Table [5J 

2.5. Line Broadening 

Observations of flare s on dMe stars 
(iHawlev fc Pettersenl 119911 : iHawlev et al.l 120031 
Kowulski et al.l 1201011 and the Sun (|Johns-Krull et alJ 
1997f) have shown that the hydrogen Balmer lines can 
become extremely broadened. These authors speculate 
the cause to be Stark broaden ing, and that concl usion 
is supported by the models of lAllred et al.l (|2006ll and 
IPaulson et al.l (1200611 . In order to test this against 
other possible sources of broadening, such as thermal or 
turbulent broadening, we have implemented a technique 
in RADYN to model Stark line broadening. The amount 
of Stark broadening depends upon the local electron 
density. Therefore, determining it can further constrain 
models of flaring stellar atmospheres. 

That the orbital angular momentum states of hydrogen 
labeled by the quantum number, 1, are not degenerate in 
the presence of an electric field perturb ation is described 
by the well-known Stark effect (e.g. iKepple fc GriemI 


119681 : IVidal etahl 119711: lSea.tonl [199011 . In stellar atmo¬ 
spheres, this perturbation is due to a net electric mi- 
croheld from fluctuations in the ambient electron, pro¬ 
ton, and ion density. The first order (linear) energy level 
shifts are proportional to the electric microfield strength 
(which has a probability distribution tY Pically modeled 
as a Holtsmark or Hooper distribution; iNavfonov et al.l 
I1999I1 . the principal quantum number, n, and a quantum 
number, q, which describes the relationship between the 
parabolic quantum numbers Oi and where Q = Q^ ~ Q -?. 
and qi + q 2 = n— | to | —1 ( Condon fc Shortl^ 1193,51 1. 
Therefore, the higher order hydrogen lines within a series 
experience a larger amount of Stark broadening. 

The general prescription for line cross-section is given 
by a Voigt profile with a damping parameter, T, that in¬ 
cludes separate terms for radiative damping. Trad, and 
for collisions among neutrals (resonance and van der 
Waals), Tres/vdw- A microturbulence of 2 km s“^ is in¬ 
cluded in the Doppler width. The Voigt profile is then 
convolved wit h the Stark p rofile to give the total line 
cross-section (jMihalasI Il978fl . However, as an approxi¬ 
mation to the more complete treatment, the Voigt pro¬ 
file damping parameter, T, can be modified to also in¬ 
clude a Stark damping term, Tg, such that T = Trad + 
Tres/vdw + Tg- This has been found to work well for 
solar hydrogen lines in the infrared (ICarlsson fc Rutt^ 
ITMai . so we have chosen this method for modeling Stark 
broadening in RADYN. 

The Stark damping parameter for hydrogen is taken 
fr om the approx imate line shape formulae in Method ^1 
of iSuttonI (|1978ll . In using this method, we are consistent 
with other NLTE radiative transfer codes such as RH 
([Uitenbroekl 1200111 . The Stark damping profile is given 
by 

rs=0.6gq{f (6) 

where gq = 0.642 for the a transition of each series and 
gq = \ otherwise, and j and i are the upper and lower 
levels of the transition, respectively. For transitions of 
helium, calcium, and magnesium, Fg = qgUe, where ga is 
a constant resulting in Stark broadening which is approx¬ 
imately three orders of magnitude less than for hydrogen. 

2.6. Boundary Conditions 

Our model flux tubes are assumed to be symmetric 
about the loop apex, so that we need only model one half 
of a full loop. The top boundary is at the loop apex where 
we have implemented a reflecting boundary condition to 
mimic incoming waves from the other side of the loop. 
The bottom boundary is below the photosphere where 
densities and pressures are very high. There we have 
implemented a simple transmitting boundary to allow 
any waves which reach that level to pass through into 
the interior. 

2.7. Initial Loops 

A focus of this work is to study how particle beams, 
which are known to be accelerated during flares, deposit 
their energy in magnetic flux loops in stellar and the 
solar atmospheres. To this end, we have generated sev¬ 
eral initial loop states with diverse lengths, temperatures, 
and densities which we will use in this study. These ex¬ 
tend from the sub-photosphere through the corona. The 
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corona is kept hot by adding a heating term, Qcor-, which 
is chosen to just balance the conductive and radiative 
losses in the upper coronal portion of the loop. Heat 
is also added to the sub-photosphere to balance radia¬ 
tive and convective flux losses there. With these heat 
sources, the loops are allowed to relax until a state of 
near equilibrium is reached. 

We have generated solar-type initial loops using three 
free parameters. These are the photospheric tempera¬ 
ture, loop-length, and coronal temperature. The coronal 
density is dependent o n loop-length and co ronal temper¬ 
ature by an RTV-type (IRosner et al.lll97§) scaling law so 
cannot be independently varied. We have chosen photo- 
spheric temperatures (i.e., the temperature in our loops 
where rsooo = 1) of 5000 K and 5800 K which correspond 
with sunspot and “quiet Sun” conditions. Loop-lengths 
were chosen to be 10 Mm for short loops and 100 Mm 
for long loops. So that we can model flaring flux tubes in 
dMe stellar atmospheres, we have also produced a model 
loop appropriate for M dwarf atmospheres. This loop has 
a higher surface gravity (562 m s“^), cooler photosphere 
(3500 K), and hotter corona (6 MK) than the solar loops. 
The parameters describing these loops are summarized in 
Table [3] and the temperature and density of the loops are 
plotted in Figure 01 

As described in Section [31 the energy deposition rate 
from particle beams depends upon the magnetic field in 
the loops. To account for this we assume a magnetic field 
strength in the photosphere of 1 kG for solar loops and 5 
kG for the M dwarf loop. Spatially-averaged observations 
of magnetic fields on active M dwarf stars have found 
that the majority of the stellar surf ace (60 — 70%) has a 
magnetic field strength of 3 — 4 kG (iSam &: LinskvllT985l : 
lSaailll99^ iJohns-Krull fc Valentilll996[l . Since these are 
spatial averages and large flares are likely to occur where 
the field is strongest, we have estimated a photospheric 
field of 5 kG in the loop. At the loop tops, the magnetic 
field strength is assumed to be 100 G for the shorter (10 
Mm) solar loops, 10 G for the longer (100 Mm) solar 
loops, and 500 G for the M dwarf loop. These values 
were chosen to be consistent with active region corona l 
magnetic field extrapolations (e.g. iTadesse et al.ll 2 ?)T^ . 
The magnetic field in all cases is assumed to exponen¬ 
tially decrease between these boundary conditions. 

3. PARTICLE BEAM HEATING 


fi = cos(a) as the pitch angle. The transit time for beam 
particles from injection source to footpoints is fast rela¬ 
tive to changes in the beam particle energy distribution 
and hydrodynamic time scales so a time-independent dis¬ 
tribution function is assumed. Thus, we will solve the 
kinetic equation for the distribution, f{E^fj,,z), where 
f{E, n, z)dEd^idz is the number density of beam par¬ 
ticles with energy between E and dE, pitch angle be¬ 
tween /r and d/i, and position between z and dz. In 
solving the kinetic equation, we include the effects of par¬ 
ticles moving at relativistic speeds. Coulomb collisions, 
synchrotron emission, pitch angle scattering, and mag¬ 
netic field gradients. However, we neglect effects due to 
plasma turbulence, electric fiel ds and self-absorption of 
synch rotron emission. We follow lMcTiernan &: Petrosian! 
(1199011 and write the Fokker-Planck equation as 

5$ d\nB d r, 2 \ 

c d r. _ 2^ ^ 

where $ = ///?, 7 = if -|-1 is the relativistic total energy 
with E measured in units of mc^, j3c is the relativistic 
ion velocity, n, B is the magnetic field strength in the 
loop, and S is a source term for particles injected at 
the loop top. C and C are coefficients which measure 
the beam energy loss rate and pitch angle diffusion from 
collisions, respectively. Similarly, S measures the loss 
rate and pitch angle diffusion due t o synchrotron emis¬ 
sion. iMcTiernan fc Petrosian! (|199Clll developed a numer¬ 
ical method to solve this equation for flare-accelerated 
electrons. We have generalized their method to solve the 
Fokker-Planck equation for both flare-accelerated elec¬ 
trons and ions. This has required altering C, C', and S 
to account for the physics which characterize ion beams. 
H ere we summariz e this theory. 

iPetrosiaiil (|1985ll calculated the energy loss and pitch 
angle scattering rates due to synchrotron radiation. 
From those results, we find the synchrotron coefficient, 
S', for both ion and electron beams to be 



Accelerated particle beams are the main source of heat¬ 
ing in the lower atmosphere during flares, so it is crucial 
that we accurately model their propagation and energy 
deposition as they move through magnetic flux loops. 
The rate of energy lost by a particle beam depends 
upon the details of the beam’s energy spectrum as well 
as the composition and ionization states of the ambi¬ 
ent plasma. We use the Fokker-Planck method coupled 
with results from our radiative hydrodynamic simula¬ 
tion described in Section to determine the distribution 
function, f{x,p,t), for flare-accelerated particles. Be¬ 
cause the magnetic field constrains the beam particles to 
move along field lines, we can replace the vector position, 
X, with the coordinate, z, which measures the distance 
along the field line. We replace the vector momentum, 
p, with the kinetic energy, E, and particle pitch angle, 
a. In the following discussion, we will simply refer to 


2 {ZefB^ 

3 {mc^)^ 


( 8 ) 


where e is the proton charge, and Z is the number of 
protons (or electrons) in the beam particle. 

C is related to the energy loss due to collisions by 


dEcoi 13 
dt c 


(9) 


IMcTiernan fc PetrosianI (1199011 assumed that the ambi¬ 
ent plasma was a “cold target” meaning that the beam 
particle’s velocity is much greater than the thermal ve¬ 
locity for all constituents of the ambient plasma. For¬ 
mally, this is expressed as Xi ^ 1 where Xi = (v/vif = 
(rni/m)E/(kTi) and rrii, Vi, and Ti are the mass, velocity 
and temperature for plasma constituent, f, and m is the 
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Table 3 
Initial Loops 


Label 

Photospheric Temperature (K) 

Loop-length (Mm) 

Coronal Electron Density (cm 

Coronal Temperature (MK) 

QS.LL.HT 

5800 

too 

3 X 10® 

3.0 

qs.ll.lt 

5800 

too 

2 X 10'^ 

1.0 

QS.SL.HT 

5800 

10 

6 X 10® 

3.0 

qs.sl.lt 

5800 

10 

6 X 10® 

1.0 

SS.LL.HT 

5000 

100 

3 X 10® 

3.0 

ss.ll.lt 

5000 

100 

2 X 10’’ 

1.0 

SS.SL.HT 

5000 

10 

6 X 10® 

3.0 

ss.sl.lt 

5000 

10 

6 X 10® 

1.0 

M dwarf 

3500 

10 

2 X 10i° 

6.0 




Figure 4. 


The temperature (a) and hydrogen number density (b) as a function of column mass for the loops described in Table [3| 


mass of the beam particle. For flare-accelerated electrons 
even at the low energy limit of ^ 10 keV, the ambient 
electron temperature would have to be greater than 100 
MK for the cold target approximation to fail. This far ex¬ 
ceeds observed temperatures for even the largest flares. 
Therefore, beam electrons can always be treated using 
cold target collision theory. With ion beams, however, 
this is not necessarily true. Flare-accelerated protons of 
^ 1 MeV (which is about the low energy detection limit 
of current instruments) have a velocity less than that of 
thermal electrons at temperature > 6 MK. Flares rou¬ 
tinely produce plasma much hotter than this. Therefore 
when modeling how ion beams interact with the ambi¬ 
ent plasma, we must employ a t heory whic h can a ccount 
for both hot and cold targets. iTrubniko^ (1196511 devel¬ 
oped this theory and found the energy loss rate of parti¬ 
cles with energy, £1, from collisions with ambient charged 
particles of species, i, to be 


dEj 

dt 


= -Evi 


( 10 ) 


where 


Vi = 2 vi 


m 


Oi- 


Erf {y/xl) - 


(l + —) 

\ TO / 


and 


VOi = 


2T:e‘^Z'^Z'f\ 
Em?c^v 


( 11 ) 

( 12 ) 


Zi is the target particle charge in units of e, and rii is the 
number density for target particle, i. Xi is the Coulomb 


logarithm, and is defined by Xi = hi{rmax/rmin)- In this 
case rmax is the mean free path length r] = v/v (|Emsliel 
119781 hereafter, referred to as E78) where v is the plasma 
oscillation frequency, is given by Tmin = h/2MiV 

where Mj is the red uced mass of the beam and target 
particle (|Huball20'Tlll . Combining these gives 


Xi = In 


Miv'^ 

h 


rrii 


^iZfe^ 


(13) 


The collisional pitch angle diffusion coefficient, C', is 
obtained from the rate that beam particles are scattered 
out of the direction of beam propagation. It is given by 


C = 


dv 

dt 


(14) 


For collisions with neutral targets (e.g., neutral hydro¬ 
gen or helium), the beam particles interact with atom¬ 
ically bound electrons. As long as the beam velocity 
is much greater than the Bohr orbital velocity, neutral 
atoms can be treated as a cold targets. Using the fine 
structure constant to approximate the orbital velocity, 
we find that the cold target approximation for collisions 
with neutral hydrogen is good for protons with energy 
> 25 keV and electrons > 0.01 keV. The energy loss rate 
for a charged beam o n a neutral target is given by (E78; 
iMott fc Massevl[l965L p. 614) 


dEn 

dt 


—27re^ ra 2 rv 

Z ZnTln^nV 

h rrie 


(15) 


where Un is the number density for the ambient neutral 
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atom, n. Zn is the number of electrons in the atom, and 
An is an effective Coulomb logarithm for collisions with 
the atom. For an atom with mean ionization energy, 

An is given by ()Evanslll95,l p. 637) 


A„, = In 


2mev'^Y 


(16) 


The total energy loss rate due to collisions is given by 

(17) 


dEcoi _ dEi dE, 


dt 


dt 


where the sums are over charged and neutral species, 
respectively. In calculating the beam energy deposition 
rate in RADYN, we consider a plasma consisting of elec¬ 
trons, protons, neutral hydrogen, neutral, singly, and 
doubly ionized helium, and singly ionized calcium and 
magnesium. These are the species most important to 
the energetics of the chromosphere where the majority of 
the beam impacts. However, since the energy loss rates 
are inversely proportional to the target particle’s mass, 
collisions with electrons -both ambient and those in neu¬ 
tral atoms- dominate the energy transfer. We use Equa¬ 
tion [10] for calculating the energy loss rate for collisions 
with charged particles and ions, and Equation [13] for col¬ 
lisions with neutral particles. We use Equations [T71 and 1^1 
to obtain C, which is needed to solve Equation |7| The 
number densities for each of these species is calculated 
in detail in the simulations. Thus, the beam energy de¬ 
position drives the simulations which alter the beam de¬ 
position in a self-consistent way. 

By solving Equation|3 we obtain the distribution func¬ 
tion, /, from which we can obtain the particle beam heat¬ 
ing rate (Qbeam) and momentum deposition rate (Abeam) 
using the following relations: 

Qbeam = (^j j t^'VEfdEdf^ (18) 

and 

Abeam = J fivaf dEdfi^ (19) 

where a = 'ymv is the relativistic momentum. Even 
though we have developed and presented here a method 
to simulate both flare-accelerated electrons and ions, in 
the remainder of this paper we will focus on the effects 
of electron beams. A detailed study of the for ion beams 
will be presented in a future paper. 


where / is the beam particle distribution function, /r is 
the particle pitch angle as described in Section [S] and aij 
is the cross section for ionization or excitation. 

Since the flare-accelerated particles have much 
larger energy than the ionization potential of hy¬ 
drogen, secondary ionizations from electrons liber- 
ated in the primary ioni zation need to be included. 
iDalgarno fc Griffmgj (1195811 have performed calculations 
to include these secondary ionizations, giving an aver¬ 
age energy of 36 eV per electron-ion pair produced in 
the complete absorption of beam particles with energy 
in the range of 200 — 1000 eV by a cold hydrogen gas. 
The amount of energy produced per electron-ion pair 
is constant for E >200 eV, and we assume that this 
extends to all electron energies above 1000 eV. Includ¬ 
ing the secondary ionizations, the non-thermal collisional 
ionization rate can be written in terms of the amount 
of energy lost by the beam to neutral hydrogen, 
obtained from Equation IT^ (iFang et al.lll993L see also 
iRicchiazzi fc Canfieldlll983r . 

^ = 36 eVxniC,(21) 

assuming that most of the neutral atoms are in the 
ground state, ni. 

A comprehensive study of no n-thermal ra tes i n 
solar flares was p r esente d by IRicchiazzi! (1198211 : 
IRicchiazzi fc Canfieldl (1198311 . who found that the H i 
rates of non-thermal collisional ionization from the 
ground state (denoted, 1-c) and collisional excitation 
from the ground to first excited state (denoted, 1-2) were 
important compared to their respective thermal rates, 
but that the non-thermal 2-c rate was negligible com¬ 
pared t o the thermal rate for the applied range of beam 
fluxes. iFang et ^ (1199311 has presented revised H i 1-c, 
1-2, 1-3, 1-4 non-therm al ionization and excitati on rates, 
which were adopted bv IKasoarova et ^ (|2009ll in their 
study of the hydrogen lines in response to electron beams. 
Notabl y, the 1-c rate is a fa c tor ^ 4.6 higher compared 
to the IRicchiazzi fc Canfieldl (1198311 1-c rate. To facili¬ 
tate comparison with th e most recent stud ies, we use the 
constants (i?^’"*) from IFang et al.l (|1993ll to derive the 
non-thermal collisional rates with hydrogen, according to 
the formula, 


/^nt _ ryH^nt 

^beam ^ 


An 


neAi + uhA, 


'Qhe 


( 22 ) 


4. PARTICLE BEAM COLLISION RATES 

During flares, the transition rates, Pij, can be very en¬ 
hanced due to direct ionizations and excitations caused 
by collisions of non-thermal beam particles with the am¬ 
bient plasma. In this section, we describe a method for 
incorporating these elevated rates in the most important 
species (i.e., hydrogen and helium) in RADYN. 

4.1. Hydrogen 

The non-thermal collisional rate from a flux of charged 
parti cles is given by the general formula (jMott fc MassevI 
119651 Chapter xvi), 

CbLm J J ^J-'VBf<7^JdEd^x ( 20 ) 


where nn is the neutral hydrogen density and = 

1.73 X 10^°, 2.94 X 10i°, 5.35 x 10®, 1.91 x 10® for the 1-c, 
1-2, 1-3, and 1-4 transitions, respectively. 

4.2. Helium 

A05 found that the ionization fraction of He ii controls 
how a flaring flux tube transitions from gentle to ex¬ 
plosive chromospheric evaporation. To more accurately 
model these dynamics, we include the non-thermal rates 
for the primary ionization of neutral helium (He i Is® —>■ 
He II Is) and s ingly ionized helium (He ii Is —>■ He ill). 
lYoungen (|1981ll calculated parameterized cross sections 
for these transitions as a function of impacting electron 
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energy. For clarity, we reproduce his result here, 

2 


uP 


1 


u 


Ail- - \ +B{1- - 


1 


u 


D 


+ C'ln(M) H-ln(M) 

u 


(23) 


where I is the ionization energy and u = E/I is the 
normalized energy of the colliding electron . A, B , C, 
and D are constants p rovid ed by lYounged (|1981ll and 
lArnaud fc Rothenflud (119851) and have values of 17.8, 
— 11.0, 7.0, —23.2 and 14.4, —5.6, 1.9, —13.3 in units 
of 10“^'^ cm^ eV^ for He i and He ii, respectively. With 
these cross sections and the beam particle distribution 
function, /, obtained from Equation [3 we calculate the 
non-thermal helium ionization rates using Equation 1201 


5. RETURN CURRENT 

Images of the footpoint locations of positively-charged 
ion beams and negatively-charged electron beams ind i- 
cate that they are not co-spatial (jHurford et al.ll200^ . 
For an electron beam without a co-streaming posi¬ 
tively charged beam of equal current density, the charge 
imbalance between the acceleration region and any 
given point along the magnetic loop leads to a macro¬ 
scopic return current electric field (iHovng et al.l 119761 : 
iKnight fc Sturroc3ll977l : Ivan den Oordll 1990(1 in a direc¬ 
tion along the beam. The return current electric field de¬ 
celerates the beam electrons while accelerating ambient 
electrontQ, which drift in the direction opposite that of 
the return current electric field towards the loop top with 
a Maxwellian distribution of speeds, ass uming that th e 
return current field is not super-Dreicer (lHolmanlll985ll . 
These drifting electrons form the return current, which 
heats the atmosphere through Ohmic (Joule) dissipation. 
Many aspects of the beam-return current-atmospheric 
system have been studie d, including pitch angle mod¬ 
ifications of the beam (|Emshd [1980) , return current 
collisional rates (IKarlickv et al.l l2004f) . return current- 
beam instabilities and subsequent particle acceleratio n 

(jKarlickv fc Kontaill2012l : lPechhacker &: Tsiklauri|l2014f) . 

and turbulent effect s on the beam-re turn current system 
(jKontar et al.l l20Tl : iXu et al.l 120131 1. The return cur- 
rent modifica tions on the classical thick target model 
(lBrownlFl971[) have been used to explain the difference 
betwe en looptop and footpoint hard X-ray s pectral in¬ 
dices (iBattaglia fc BenjllMOSl : IXu et al.ll2013ll . 

Recently. iHolmanI (l2012f) "derived return current heat¬ 
ing rates in flaring coronal conditions. After the first 
short moments of beam pr opagation, a relative ly steady- 
state is quickly attained ([van den OordlfTfi^ in which 
the magnitude of the return current density is equal the 
beam current density. The return current electric field 
can be determined from Ohm’s law: Ere = fl^heam, where 
rj is the classical Spitze r resistiv i ty and Jbeam is the re¬ 
turn current density. IHolmanI (j2012l) determined the 
return-current plasma heating rate for a given electron 
beam flux self-consistently accounting for Joule heat¬ 
ing and the reduction of beam electrons due to (non- 
collisional) thermalization caused by the return current 


^ Ambient protons are also accelerated, but the drift velocity is 
much lower due to their larger mass. 


electric field. However, this treatment does not self- 
consistently account for the change in flux of beam elec¬ 
trons due to Coulomb collisions with the ambient plasma 
and is therefore an overestimate, especially in the chro¬ 
mosphere where Coulomb collisions are large. In the 
corona, however, there are relatively few collisions, and 
this approximation works quite well. Additionally, re¬ 
turn current heating is likely largest in the corona, since 
it is there that beam current densities are highest. There¬ 
fore, to account for return current heating in fl ares, we 
have chosen to incorporate the formalism of IHolmanI 
(|20I2f) into RADYN. The result for the volumetric re¬ 
turn current heating rate from a power-law electron beam 
spectrum with power-law index, S, and cutoff energy, Ec, 
is given by 




rie^Ee^ ^(5 -I- 

.X + 


1-2(5 


X < Xrc 


X > Xrc 


(24) 


where x is distance from the loop top and Ec is the 
beam particle flux given by / / ^vfdEd^. Beam elec¬ 
trons are assumed to be thermalized and removed from 
the beam when their energy is Etherm = 2.5kT. Xrc 
is the distance at which the lowest energy electrons in 
the beam are thermalized, an d V(x) is the e lectric po¬ 
tential energy as described in IHolmanI (I20I2D . A more 
complete treatment which includes return current losses 
directly in the Fok ker-Planck equation (Eqn. O, simi¬ 
lar to work done by [Z har kova fc Gordovskvvi (120051 1 and 
Dobranskis fc Zharkovar(l20I4ll . is outside the scope of 
this paper but will be presented in a future paper. 


6 . PARAMETER STUDY 

The energy distribution of flare-accelerated particles 
can vary greatly between particular flares. Since the im¬ 
pact location strongly depends on the particle distribu¬ 
tion, that distribution is very important in determining 
how stellar atmospheres respond to flare heating. To ex¬ 
plore the range of possibilities, we have performed a pa¬ 
rameter study varying the injected electron beam spec¬ 
trum, pitch-angle distribution, and the initial loop con¬ 
ditions into which these particles impact. 

Flare-accelerated particles are typically observed to 
have a power-law energy distribution with a cutoff en¬ 
ergy, Ec, below which relatively few particles are accel¬ 
erated. Cutoff energies as high as 250 keV have been 
observed (R. Schwartz, personal communication, 2015). 
But in many flares, the cutoff energy cannot be directly 
determined, since thermal X-ray emission in the range 
below 10 — 20 keV from the flare-heated plasma over¬ 
whelms the non-thermal bremsstrahlung emission. In 
these flares, it is only possible to determine an upper 
limit to the cutoff energy. To explore how the atmosphere 
responds to beams of particles with energies below this 
upper limit, we have chosen a lower limit of 5 keV. There¬ 
fore, to account for the full range in parameter space, we 
have varied the cutoff energy between 5 — 500 keV. The 
slope of the power-law is given by the spectral index, S. 
Typical power law indices range from 3 — 9, but in this 
study we have varied it between 2.5 — 10 to ensure that 
we have spanned the full range of possible values. We 
have also varied the pitch angle distribution with which 
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Table 4 

Parameter Study 


Parameter 

Range 

Ec 

5 - 500 keV 

5 

2.5 - 10 

^0 

0.1 — 1.0, isotropic, beamed 

Loop Conditions 

L’hose listed in 'l‘able[H| 


the flare-accelerated particles were injected. We modeled 
this distribution as a Gaussian centered around the flux 
loop axis. Thus, it has the form, , where p,o 

is a parameter which controls the width of the Gaussian 
distribution. In addition to the Gaussian distribution, we 
have included simulations in which the particles are fully 
beamed, i.e., they all start with a 0° pitch-angle, and 
simulations with isotropic distributions. We have mod¬ 
eled the effects of beam impact varying these parameters 
on each of the loops in Table |3l The parameters are sum¬ 
marized in Table 01 To perform this study, we ran more 
than 11,000 simulations. 

6 .1. Penetration Depths in Flux Loops 

First, we consider how varying Ec and S affect the lo¬ 
cation of beam impact. Figure [S] shows a few represen¬ 
tative examples of the heating rate due to collisions with 
beam particles as a function of hydrogen column den¬ 
sity in the QS.SL.HT loop. Glearly the heating extends 
over a range of column densities. We have determined an 
average value indicated by the dashed lines. In the fol¬ 
lowing discussion, the term penetration depth (in units 
of column density) refers to this average. Figured shows 
the beam penetration depth as a function of cutoff en¬ 
ergy and the power-law index in the QS.SL.HT and M 
dwarf flux tubes. A striking result of this analysis is 
how strongly dependent penetration depth is on S. Low 
values of 6 imply a significant number of higher energy 
particles even for spectra with low cutoff energies. For 
example, an electron distribution with Ec = 10 keV and 
S = 3.0 penetrates as deeply as one with Ec = 70 keV 
and 6 = 5.0. Interestingly, for distributions that have 
both very low values of Ec (< 8 keV) and high 6 (> 6) 
much of the energy is deposited directly in the corona. 

Many previous studies fe.g.. lAbbett fc HawlevI 119991 : 
I Allred et al.ll2005l . 120061: iGheng et al.ll2010h implemented 
the analytic expression for beam penetration depth de¬ 
rived b y E78. That expression includes non-uniform ion¬ 
ization (|Hawlev fc Fisheill 19941) but does not include rel¬ 
ativistic effects. Since the E78 expression has been so 
widely used, it is informative to compare it with the more 
complete treatment that we described in Section [31 In 
Figure [71 we plot the ratio of beam penetration depths 
obtained from E78 to that derived by solving the Fokker- 
Planck equation. We find E78 works quite well in the the 
range Ec < 40 keV and <5 > 4.5 but becomes increasingly 
inaccurate for higher Ec and lower S. In this regime, E78 
predicts a penetration depth as much as 7 times greater. 
As the cutoff energy increases, relativistic effects become 
significant. Due to the Lorentz length contraction, rel¬ 
ativistic electrons experience a higher ambient density 
than would be expected from classical theory resulting 
in less penetration. Relativistic effects become less pro¬ 
nounced with increasing J, since there are fewer high 


Height (Mm) 

9.40 1.20 1.11 0.90 0.65 0.39 0.12 



Figure 5. The beam heating rate for several different parameters 
in the QS.SL.HT loop. In the top panel, /io is varied while holding 
Ec = 20 keV and (5 = 5. In the middle panel, Ec is varied while 
holding (5 = 5 and fiQ = 0.1. In the bottom panel, 5 is varied while 
holding Ec = 20 keV and {iq = 0.1. In each panel, the dotted line 
indicates the temperature structure in the QS.SL.HT loop. 

energy electrons in the distribution. 

Figure [8] shows the ratio of the penetration depth for 
highly beamed non-thermal electrons to an isotropic in 
pitch angle distribution. As expected, the beamed elec¬ 
trons penetrate more deeply. As 6 increases this effect 
becomes more pronounced. This is because higher en¬ 
ergy particles undergo more collisions before stopping, 
so their initial pitch angle is less important in determin¬ 
ing their final stopping depth. A larger S means fewer 
high energy particles are in the distribution, so the initial 
pitch angle has a more pronounced effect. 

Figure |3| compares the beam penetration depths as a 
function of cutoff energy and spectral index for several 
loops. The left panel shows the ratio of the penetration 
depths for beams propagating in QS.SL.HT to QS.SL.LT. 
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Average Penetration Depth 



10 100 
Cutoff Energy (keV) 


T(kK) 



4.6- -21.8g 


4.7- -21.0 I 
Q 



Average Penetration Depth 



10 100 
Cutoff Energy (keV) 


T{kK) 

3.0^H23.5 

3.iH^2-7 

4.2 

5.6 

6.5 

3.4 MK^H 19.3 
5.7 MK^I 18.5 


21.8 i 
21.0 i 

Q 

c 

20.2-i 


Figure 6. The penetration depth as a function of cutoff energy and spectral index, S, in the QS.SL.HT atmosphere (left) and M dvt^arf 
atmosphere (right). 


Ratio of E78 to FP Penetration Depth 



10 100 
Cutoff Energy (keV) 


Figure 7. Ratio of the beam penetration depth calculated using 
our Fokker-Planck technique to that of E78 in the QS.SL.HT loop. 


Ratio of Beamed to Isotropic Penetration Depth 



10 100 
Cutoff Energy (keV) 


Figure 8. Ratio of the penetration depth for a highly beamed 
to isotropic injected distribution. Calculated using the QS.SL.HT 
initial atmosphere. 

For low values of the cutoff energy, the beam penetrates 
to a column mass approximately twice as great in the 
QS.SL.lt loop. In this loop, the coronal temperature 
and density are lower than in QS.SL.HT. It is the low 


energy electrons that are most affected by this difference. 
The middle panel compares beam penetration depths in 
QS.LL.lt with QS.SL.lt. In this case, beams with low- 
cutoff energy penetrate less deeply in QS.SL.LT. Since 
QS.SL.lt is a shorter loop than QS.LL.LT, the coronal 
density is greater resulting in less penetration for the low¬ 
est energy electrons. The right panel compares QS.SL.LT 
with SS.SL.lt. These loops are nearly identical in their 
coronal and upper chromospheric regions. They differ 
significantly in the photosphere, but very few electrons 
are able to penetrate to that depth so beam penetra¬ 
tion depths are similar in these loops in the range of 
cutoff energies and spectral indices studied here. How¬ 
ever, since they have very different photospheric temper¬ 
atures, these loops predict significantly different levels of 
white light emission. This difference will be important in 
studying the white light emission predicted from flaring 
loops. 

6 .2. Collision Rates 

During flares, direct excitations and ionizations by 
beam particles can be much larger than thermal rates in 
the region of beam impaclO. To understand the relative 
importance of collisional ionizations by flare-accelerated 
particles, it is informative to compare them with the ther¬ 
mal rates. In Figure [TUI plot the ionization rates, 
niC^^am^ in several loops. Here is the ground state 
number density. We find the ionization rates increase 
dramatically around beam impact. For example, in the 
solar loops the H i 1-c rates peak at 2 x 10^^ s”^. To put 
this into perspective, the thermal rates in this region are 
- 10^ s-i so some 10 orders of magnitude smaller. These 
results are relatively independent of the initial loop con¬ 
ditions, with the exception of the M dwarf loop. 
increases deeper in the M dwarf loop compared to the so¬ 
lar loops because in the former the chromosphere forms 
deeper. 

It is useful to model how these collision rates vary 
with injected electron beam spectra. Figure ITT] shows the 

^ Even though these excitations are important, throughout the 
remainder of this section, we will focus exclusively on the ioniza¬ 
tions. This is because the excitation rates can be obtained simply 
from scaling the ionization rates by the Rtr.nt coefficients as seen 
in Equation 1221 
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Figure 9. Ratio of the beam penetration depth in the QS.SL.HT (left), QS.LL.LT (middle), and SS.SL.LT (right) atmospheres to the 
penetration depth in the QS.SL.LT simulation as a function of cutoff energy and spectral index. 



Figure 10. The beam induced collisional rates, riiCg^, for H 1 1-c 
(solid line), He I 1-c (dashed line), and He ll 1-c (dotted line) in 
the QS.SL.HT (black), QS.SL.LT (red), QS.LL.HT (green), and M 
dwarf (blue) loops. In each case the injected electron beam is for 
Ec = 20 keV and <5 = 5. The He ll 1-c rates have been scaled by a 
factor of 1000 so that they will fit in the plot. 

peak non-thermal collisional ionziation rates for electron 
beams injected into the QS.SL.HT loop as a function of 
Ec and 6. H i 1-c and He i 1-c rates have similar de¬ 
pendence on these parameters, peaking at lower Ec and 
higher S. Qbeam is narrower and higher peaked in this 
regime (see the middle and bottom panels of Figure [5]) 
resulting in a higher peak collision rate. However, the 
collision rates drop quickly for very low Ec- There much 
of the beam is stopped in the transition region and corona 
where rii is low. The He ii 1-c rates rapidly decrease for 
increasing Ec- A beam with larger Ec penetrates to a 
deeper, cooler, and hence lower He ii density region of 
the atmosphere. 

6.3. Return Current Heating in the Corona 

To understand how return currents are likely to affect 
flare dynamics, it useful to compare their heating rates 
{Qrc) to the heating rates produced directly by collisions 
(Qbeam)- In the following discussion, it should be kept in 
mind that we are comparing only coronal heating rates - 
not those in the chromosphere where Qbeam dominates. 
A comparison of these rates in several loops is shown in 
Figure [121 In the QS.SL.HT loop, Qrc and Qbeam are 
similar in size. However, in all other cases, Qrc domi¬ 


nates by at least an order of magnitude, being stronger 
in the cooler loops. This is because the resistivity is 
strongly and inversely dependent on the coronal temper¬ 
ature. Hot loops have much less resistance so are heated 
by Qrc less. Therefore, we speculate that return cur¬ 
rents are likely to be most important in the early phases 
of flares before the corona has been heated to very high 
temperatures. This figure also compares Qrc for the case 
of highly beamed electrons (solid lines) and isotropically 
injected electrons (dotted lines). These rates are similar 
in all cases indicating that the pitch-angle distribution of 
the injected electrons makes little difference in Qrc- 

From the results of this parameter study, we can also 
compare how Qrc varies with Ec and <5. Figure [13] shows 
the volumeteric heating rate due to return current in the 
corona of the QS.SL.HT loop as a function of these pa¬ 
rameters. The same energy flux (10^° erg cm“^ s“^) was 
injected for every simulation in this study, so simulations 
with higher Ec have fewer total injected electrons. Re¬ 
turn current heating scales with the particle flux (Ec) so 
fewer particles results in less heating. The energy per 
particle factors into Ec only through its velocity, so the 
return current heating rate is only weakly dependent on 
( 5 . 

7. CONCLUSIONS 

We have developed a computational framework capa¬ 
ble of modeling the evolution of magnetic flux loops in 
response to flare heating. This model can be applied to 
loops on the Sun and on dMe flare stars, and it includes 
many of the physical processes most important to flare 
dynamics. Our code solves the optically-thick, non-LTE 
radiative transfer equation for many transitions impor¬ 
tant in the chromosphere where much of the flare energy 
is deposited allowing us to model flare emission which 
can be compared directly to observations. We have im¬ 
plemented a method for solving the Fokker-Planck equa¬ 
tion describing how flare-accelerated particles interact 
with the ambient plasma and have used that to model 
the heating, momentum deposition, collisional excita¬ 
tions and ionizations, and return currents produced by 
these particles. 

With our model, we have performed a parameter study 
to understand how these phenomena depend on the en¬ 
ergy spectrum of a beam of flare-accelerated particles. 
We have assumed a power-law energy distribution for 
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Figure 11. Collisional ionziation rates for H i (top), He i (middle), 
and He ii (bottom) produced by non-thermal particles injected at 
the top of the QS.SL.HT loop. 

these particles and found that their penetration depth 
is strongly dependent on the beam cutoff energy and 
power-law index and more weakly dependent on their 
initial pitch-angle distribution. For distributions with 
low power-law index or high cutoff energy, relativistic ef¬ 
fects become important causing a penetration to depths 
as much as 7 times less than predicted by classical the- 
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Figure 12. Volumetric heating rates for the return current 
(solid and dotted lines) and electron beam (dashed lines) in the 
QS.SL.HT (black), QS.SL.LT (red), QS.LL.HT (green), QS.LL.LT 
(orange), and M dwarf loops (blue). The solid and dashed lines 
indicate return current heating for beamed and isotropic distri¬ 
butions of electrons, respectively. These rates are all calculated 
assuming Ec = 20 keV and 5 = 5. 
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Figure 13. Volumetric return current heating rate as a function 
of Ec and 6 in the QS.SL.HT loop. 


ory. Varying the conditions of the loops (i.e., the coronal 
temperature and density and photospheric temperature) 
has a small effect on the penetration depth for distribu¬ 
tions with low cutoff energy. The collisional ionizations 
produced by impacts from these flare-accelerated parti¬ 
cles also are strongly dependent on cutoff energy and 
power-law index. They can be 10 orders of magnitude 
greater than thermal collision rates. The beam particles 
also induce a return current which deposits heat in the 
corona. We have found that, depending on coronal con¬ 
ditions, return current heating in the corona can be more 
than an order of magnitude greater than coronal beam 
heating. 

In this parameter study, we have focused on under¬ 
standing the effects of flare-accelerated electrons. In fu¬ 
ture work, we will extend this to flare-accelerated ions 
which are also known to be important to flare energetics. 
This work has focused on the initial response of flaring 
loops. However, a key feature of our framework is its abil¬ 
ity to model the dynamical evolution of these loops. In 
future work, we will present a parameter study exploring 
their evolution in response to flare-accelerated particles. 
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